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Abstract 

We propose a method to study lattice QCD at finite temperature (T) 
and chemical potential (/i). We compare the method with direct results 
and with the Glasgow method by using nf=A QCD at Im(;u)7^0. We 
locate the critical endpoint (E) of QCD on the Re(/x)-T plane. We use 
n/=2+l dynamical staggered quarks with semi-realistic masses on L(=4 
lattices. Our results are based on 0(10^ — 10^) configurations. 

QCD at finite T and/or /i is of fundamental importance, since it describes 
relevant features of particle physics in the early universe, in neutron stars and 
in heavy ion collisions. Extensive experimental work has been done with heavy 
ion collisions at CERN and Brookhaven to explore the /z-T phase boundary (cf. 
[H]). Note, that past, present and future heavy ion experiments with always 
higher and higher energies produce states closer and closer to the T axis of the 
/i-T diagram. It is a long-standing question, whether a critical point exists on 
the ii-T plane, and particularly how to tell its location theoretically 0. 

Let us discuss first the /i=0 case. Universal arguments and lattice simu- 
lations 1^ indicate that in a hypothetical QCD with a strange (s) quark mass 
{rris) as small as the up (u) and down (d) quark masses {mu,d) there would be a 
first order finite T phase transition. The other extreme case (n/=2) with light 
u/d quarks but with an infinitely large there would be no phase transition 
only a crossover. Observables change rapidly during a crossover, but no singu- 
larities appear. Between the two extremes there is a critical strange mass (m^) 
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at which one has a second order finite T phase transition. Staggered lattice 
results on Lt=4: lattices with two light quarks and around the transition T 
{uf =2+1) indicated that is about half of the physical m^. Thus, in the 
real world we probably have a crossover. 

Returning to a non-vanishing /i, one realizes that arguments based on a 
variety of models (see e.g. 0, |^) predict a first order finite T phase transition 
at large /i. Combining the fi = and large /i informations an interesting 
picture emerges on the /i-T plane. For the physical rris the first order phase 
transitions at large fi should be connected with the crossover on the /i = 
axis. This suggests that the phase diagram features a critical endpoint E 
(with chemical potential fiE and temperature Te), at which the line of first 
order phase transitions {ji> he and T < Te) ends 0. At this point the phase 
transition is of second order and long wavelength fiuctuations appear, which 
results in (see e.g. |^) consequences, similar to critical opalescence. Passing 
close enough to {(ieiTe) one expects simultaneous appearance of signatures 
which exhibit nonmonotonic dependence on the control parameters [^, since 
one can miss the critical point on either of two sides. 

The location of E is an unambiguous, non-perturbative prediction of QCD. 
No ah initio, lattice analysis based on QCD was done to locate the endpoint. 
Crude models with = cx) were used (e.g. 0) suggesting that /i^ ~ 700 MeV, 
which should be smaller for finite m^. The goal of our exploratory work is to 
show how to locate the endpoint by a lattice QCD calculation. We use full 
QCD with dynamical nf =2+1 staggered quarks. 

QCD at finite /i can be given on the lattice however, standard Monte- 
Carlo techniques fail. At R.e{fi)^0 the determinant of the Euclidean Dirac 
operator is complex, which spoils any importance sampling method. 

Several suggestions were studied in detail to solve the problem. We list a 
few of them (for a recent review see Ref. ^^)■ 



In the large gauge coupling limit a monomer-dimer algorithm was used [12 



For small gauge coupling an attractive approach is the "Glasgow method" |TB 



in which the partition function is expanded in powers of exp(yu/T) by using 
an ensemble of configurations weighted by the /i=0 action. After collecting 
more than 20 million configurations only unphysical results were obtained: 
a premature onset transition. The reason is that the /i=0 ensemble does not 
overlap sufficiently with the states of interest. Another possibility is to separate 
the absolute value and the phase of the fermionic determinant and use the 



former to generate configurations and the latter in observables [14 



At imaginary n the measure remains positive and standard Monte Carlo 
techniques apply. The canonical partition function can be obtained by a 



Fourier transform O, IT6|. In this technique the dominant source of errors 



is the Fourier transform rather than the poor overlap. One can also use the 
fact that the partition function away from the transition line should be an 
analytic function of n, and the fit for imaginary fi values could be analytically 



continued to real values of fi At T sufficiently above the transition, both 
real and imaginary fi can be studied by dimensionally reducing QCD ||18 
Hamiltonian formulation may also help studying the problem [|19|] . 



new method 
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Fig. 1: Schematic difference between the present and the Glasgow methods. 

An attractive approach to alleviate the problem is the "Glasgow method" 
(see e.g. Ref. |jl3|) in which the partition function [Z) is expanded in powers of 
exp(/i/T) by using an ensemble of configurations weighted by the /i=0 action. 
After collecting more than 20 million configurations only unphysical results 
were obtained. The reason is that the /i=0 ensemble does not overlap enough 
with the finite density states of interest ||2^ . 

We propose a method to reduce the overlap problem and determine the 
phase diagram in the /i-T plane (for details see |21[|). The idea is to produce 
an ensemble of QCD configurations at /x=0 and at T^. Then we determine 
the Boltzmann weights of these configurations at /i 7^ and at T low- 
ered to the transition temperatures at this non-vanishing fi. Since transition 
configurations are reweighted to transition ones a much better overlap can be 
observed than by reweighting pure hadronic configurations to transition ones 
[0 . Since the original ensemble is collected at /i=0 we do not expect to be able 
to decsribe the physics of the large /i region with ie.g. exotic colour supercon- 
ductivity. Fortunately, the typical fi values at present heavy ion accelerators 
are smaller than the region we cover (for the freeze-out condition see e.g. p^). 

After illustrating the applicability of the method we locate the critical 
point of QCD. (Multi-dimensional reweighting was successful for determining 



the endpoint of the hot electroweak plasma |22| e.g. on 4D lattices.) 



Let us study a generic system of fermions ip and bosons 0, where the fermion 
Lagrange density is ipM{(j))ilj. Integrating over the Grassmann fields we get: 



Z(a) 



V(f)exp[—Shos{ci, 0)] det M(0, a), 



(1) 



where a denotes a set of parameters of the Lagrangian. In the case of staggered 
QCD a consists of /3, and fi. For some choice of the parameters a=ao 
importance sampling can be done (e.g. for Re(/i)=0). Rewriting eq. (|l|) 



(2) 



Z(a) = I r'0exp[-S'6os(«o, 0)] det M(0, ao) 

detM(0,Q!) 1 



exp[-Sbosia, 0) + Sbosiao, 0)] 



detM(0,ao) 



We treat the curly bracket as an observable (measured on each configuration) 
and the rest as the measure. Changing only one parameter of the ensemble 
generated at provides an accurate value for some observables only for high 
statistics. This is ensured by rare fiuctuations as the mismatched measure 
occasionally sampled the regions where the integrand is large. This is the 
overlap problem. Having several parameters the set ao can be adjusted to get 
a better overlap than obtained by varying only one parameter. 
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Fig. 2: 'tpip as a function of Im(/i), for direct results (squares), our technique (crosses) and 
Glasgow-type reweighting (dots). 

The basic idea of the method as applied to dynamical QCD can be summa- 
rized as follows. We study the system at Re(/i)=0 around its transition point. 
Using a Glasgow-type technique we calculate the determinants for each con- 
figuration for a set of /i, which, similarly to the Ferrenberg-Swendsen method 



W[ , can be used for reweighting. The average plaquette values can be used to 



perform an additional reweighting in (3. Since transition configurations were 
reweighted to transition ones a much better overlap can be observed than by 
reweighting pure hadronic configurations to transition ones as done by the 
Glasgow-type techniques. The differences between the two methods are shown 



in Figure 1. Moving along the transition line was also suggested by Ref. |T6|. 
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Fig. 3: Im(/3(J°) as a function of the chemical potential. 

We have directly tested these ideas in ri/=4 QCD with 777.5=0. 05 dynamical 
staggered quarks. We first collected 1200 independent V=4-6^ configurations 
at Re(/i)=lm(/i)=0 and some [3 values and used the Glasgow-reweighting and 
also our technique to study Re(/i)=0, Im(/i)7^0. At Re(/i)=0, Im(/i)7^0 direct 
simulations are possible. After performing these direct simulations as well, a 
clear comparison can be done. Figure 2 shows the predictions of the three 
methods for the average quark condensates at /3=5.085 as a function of Im(/i). 
The predictions of our method agree with the direct results, whereas for larger 
Im(//) the predictions of the Glasgow method are by several standard devia- 
tions off. We expect that our method can be applied at Re(/i)7^0. 

In QCD with 77^ staggered quarks one should change the determinants to 
their 77/ /4 power in our two equations. Importance sampling works also in this 
case at some [3 and at Re(/i)=0. Since det M is complex an additional problem 
arises, one should choose among the possible Riemann-sheets of the fractional 
power in eq. (0). This can be done by using the fact that at fi=fiw the 
ratio of the determinants is 1 and it should be a continuous function of /i. 

In the following we keep /i real and look for the zeros of Z for complex f3. 
At a first order phase transition the free energy oc log Z{f]) is non-analytic. 
A phase transition appears only in the V— 00 limit, but not in a finite V. 
Nevertheless, Z has zeros at finite V, generating the non-analyticity of the free 
energy, the Lee- Yang zeros These are at complex parameters (e.g. 



For a system with first order transition these zeros approach the real axis as 
V-^ oo by a 1/V scahng. This V— > oo hmit generates the non-analyticity of 
the free energy. For a system with crossover Z is analytic, and the zeros do 
not approach the real axis as V— > oo. 
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Fig. 4: The T-/i diagram. Direct results are given with errorbars. Dotted hne shows the 
crossover, sohd hne the first order transition. The box gives the uncertainties of the endpoint. 

At Tt^O we used Lt=4, Ls=4,6,8 lattices. T=0 runs were done on 10^- 
16 lattices. mu,rf=0.025 and ms=0.2 were our bare quark masses. At T 7^ 
we determined the complex valued Lee- Yang zeros, /3o, for different V-s as 
a function of ji. Their V— > 00 limit was given by a (SoiV) — + C/^ 
extrapolation. We used 14000, 3600 and 840 configurations on Ls=4,6 and 8 
lattices, respectively. Im(/?^) is shown on Figure 3 as a function of fi. For 
small fi-s the extrapolated lm{f3^) is inconsistent with a vanishing value, and 
predicts a crossover. Increasing the value of Im(/5^) decreases, thus the 
transition becomes consistent with a first order phase transition (overshooting 
is a finite V effect). Our primary result is /lend — 

0.375(20). 

To set the physical scale wc used a weighted average of Rq, rrip and ^/a. 
Note, that (including systematics due to finite V) we have {Rq -171^^) = 0.73(6), 
which is at least twice, rriu^d is at least four times as large as the physical values. 

Figure 4 shows the phase diagram in physical units, thus T as a function 
of /ib, the baryonic chemical potential (which is three times larger then the 
quark chemical potential). The endpoint is at Te = 160 ± 3.5 MeV, — 
725 ± 35 MeV. At hb=0 we obtained = 172 ± 3 MeV. 

We proposed a method -an overlap improving multi-parameter reweighting 
technique- to numerically study non-zero /i and determine the phase diagram 
in the T-ji plane. Our method is applicable to any number of Wilson or 
staggered quarks. As a direct test we showed that for Im(/.()7^0 the predictions 



of our method are in complete agreement with the direct simulations, whereas 
the Glasgow method suffers from the well-known overlap problem. We studied 
the n-T phase diagram of QCD with dynamical nf=2+l quarks. Using our 
method we obtained Te;=160±3.5 MeV and yU£;=725±35 MeV for the endpoint. 
Though /is is too large to be studied at RHIC or LHC, the endpoint would 
probably move closer to the fi=0 axis when the quark masses get reduced. At 
/i=0 we obtained Tc=172±3 MeV. More work is needed to get the final values 
by extrapolating in the R-algorithm and to the thermodynamic, chiral and 
continuum limits. The details of the presented results can be found in [pTj . 
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